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THE  FACTORIZATION  METHOD  FOR  TWO  POINTS  BOUNDARY  VALUE  PROBLEMS 
FOR  ODE'S  AND  ITS  RELATION  TO  THE  FINITE  DIFFERENCE  METHOD 


I.  Babuska^*  and  V.  Majer^ 


1.  INTRODUCTION 

Finite  difference  and  finite  element  methods  for  solving  two  point 
boundary  value  problems  for  systems  of  ordinary  differential  equations  consist 
of 

a)  a  discretization  procedure  which  transforms  the  original  problem 
into  a  family  of  finite  dimensional  systems  of  algebraic  equations 
parametrized  by  the  mesh  size  h,  and 

b)  a  solution  procedure  for  the  systems  of  algebraic  equations. 

For  linear  boundary  value  problems  the  algebraic  equations  are  linear  and  step 
b)  reduces  to  the  selection  of  a  matrix  reduction  scheme.  In  this  paper  we 
consider  only  direct  (elimination)  methods  of  matrix  reduction. 

By  these  two  steps,  taken  together,  the  original  problem  is  transformed 
into  4  sequential  numerical  process  (§5)  which  depends  on  the  mesh  parameter 
h.  A  comptetft  analysis  of  the  numerical  procedure  must  consider  this 
underlying  numerical  process,  not  merely  the  discretization  etep  a).  In  this 
paper  m  carry  out  su*h  a  complete  analysis  for  a  model  singular  perturbation 
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problem  of  turning  point  type  (§2)  studied  by  H.  0.  Kreiss  et  al.  in  [1],  We 
show  (§4)  that  the  numerical  process  converges,  as  h  0,  to  the  solution  of 
initial  value  problems  for  certain  differential  equations.  These  limiting 
equations  are  the  closure  [2]  of  the  process.  Thus  it  is  possible  to 
interpret  the  numerical  process  as  a  special  method  (of  low  order)  for  solving 
these  initial  value  problems. 

This  fact  suggests  that  one  should  study  directly  a  class  of 
transformations  of  the  original  boundary  value  problem  into  systems  of  initial 
value  problems,  called  factorizations,  which  class  includes  the  above- 
mentioned  closure.  We  address  this  question  in §3,  where  we  single  out  those 
transformations  of  this  class  which  are  stable  in  a  precisely  defined  sense, 
and  which  therefore  can  be  solved  by  proper  numerical  methods.  One  such 
method,  of  course,  would  lead  to  the  identical  numerical  process  as  that 
stemming  from  the  finite  element  or  finite  difference  discretizations. 

Although  we  restrict  ourselves  here  only  to  a  linear  model  problem,  the 
results  described  here  hold  for  general  systems  of  boundary  value  problems  and 
can  be  applied  also  to  nonlinear  problems  [3], 
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2.  MODEL  PROBLEM 


2.1.  Statement  of  the  problem 


Consider  the  scalar  second’-order  boundary  value  problem 


(2.1*0 


(2  .m 


d 

— 2  U<'8^  *  3?  (a(s)u(s))  +  b(s)u(s)  +  g(s),  s  €  [o,d*] 

ds 


v(a)  *  y. 


«#(€r*)  »  y* .. 


The  functions  a,  t,  %  &te  smooth  osn  s*»d  0  <  e  «  1  is  a 

constant  parameter,  ISsjk  solsHSiana.  of  (2.1)  can  exhibit  boundary  or  interior 
layers  in  whi<$*  0%*  solution  is  rapidly  changing,  A  discussion  of  the  range 
of  possible  behaviors  may  be  found  in  [1], 

We  seek  the  values  O-f  an  approximate  solution  of  (2.-1)  on  a  partition 


x:  0  *  *n  <  ffj  <  •••  <  °m  * 


of  target  points  of  [o,o*].  We  -s'ieo  like  tc*  determine  the  location  of 

tfui  interior  or  boundary  layers  and  perhaps  to  resolve  them. 

2.2.  Reformulation  as  a  First  Order  System 

It  is  possible  to  recast  (2.1)  as  a  first  order  system  in  a  number  of 
ways.  We  consider  one  of  them  here.  Let  v'  «  bu  +  g.  Then 


u'(s) 


^u(s)  +  ±v(s), 


(2.2a) 


(2.2b) 


v'(s)  -  b(s)u(s)  +g(s). 


u(a)  ■  y,  v(o*)  -  y*> 


1 


in  which  we  have  used  —  ■  ' .  Written  in  matrix  form  (2.2)  becomes 


’u'(s)' 

*a(s)/e  1/c- 

u(  s)‘ 

■  0  - 

(2.3a) 

m 

+ 

_v'(v)_ 

_b(  s)  0_ 

v(s)_ 

_g(  s)_ 

subject  to 


"  u(a)  ' 

■u(o*r 

(2.3b) 

[1  0] 

-  Y. 

[1  0] 

. v(a)  . 

_v(  a* ) 

Equation  (2.3)  is  but  a  special  case  of 


'u'(s)- 

-a(s)  c(s)‘ 

'  u(s)‘ 

■f(s)- 

(2.4a) 

- 

+ 

,v'(s)_ 

_b(s)  d(s)_ 

_  v(s). 

_g(  s). 

subject  to 


"u(a)‘ 

"u(o*)“ 

(2.4b) 

[a  6] 

=  V 

[a*  3*] 

.v(a)_ 

_v(a*). 

The  form  (2.4)  encompasses  all  the  various  reductions  of  (2.1)  to  first  order 
form.  In  the  sequel,  it  is  assumed  that  the  solution  of  (2.4)  exists. 
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3.  INDIRECT  SOLUTION  OF  THE  MODEL  PROBLEM 

The  indirect  method  of  solution  of  (2.4)  is  based  upon  the  propagation 
across  the  interval  [o,o*]  of  the  boundary  conditions  (2.4b)  in  a  manner 
consistent  with  the  original  equation.  For  this  we  require  certain  auxiliary 
initial  value  problems. 


3.1.  Auxiliary  Initial  Value  Problems 

Consider  the  initial  value  problems 


(3.1) 


(3.2) 


(3.3) 


(3.4) 


(3.5) 


<j>'(s)  =  -a(s)<|>(s)  -  b(s)i)»(s)  +  z(s)ij)(s), 


$(o)  *  a. 


ijr(s)  =  -c(s)(f)(s)  -  d(sH(s)  +  z(s)ip(s), 


♦(c)  =  8, 


/'(s)  -  f(s)<(>(z)  +  g(s)i{>(s)  +z(s)x(s). 


x(o)  -  y » 


♦*'(s)  “  -a(s)$*(s)  -  b(s)ty*(s)  +  z*( s) 4»*( s) , 


>(o*)  -  o*, 


t|>*'(s)  -  -c(s)$*(s)  -  b( s)\|;*( s)  +  z*(s)<J/*(s) , 


♦*(o)  -  B*. 
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X*(s)  =  f  ( s>  <)>*(  s)  +  g(s)\|>*(s)  +  z*(s)x*(s), 

(3.6) 

x*(a*)  =>  y* , 

in  which  the  conditioning  functions  s  ♦  z(s),  z*(s)  are 

the  moment,  otherwise  arbitrary.  Problems  (3.1)-(3.3)  are 
while  (3.4 )— ( 3.6)  are  solved  backward. 

Theorem  3.1.  Suppose  that  (3.1 )— ( 3.3)  have  been  solved  on  (o,5)  and  (3.4)- 
(3.6)  have  been  solved  on  Then  for  s  €  [cr,£]  fl  (5*, a*]  ,  and 

[u(s)v(s)]T  satisfying  (2.4), 


>(s) 

ty(s) 

"u(s)“ 

'x(s)  ■ 

_4>*(s) 

t*(s). 

.v(s)_ 

-X*(s)_ 

Corollary  3.2.  If  s  >  [u(s)  v(s)]T  is  the  unique  solution  of  (2.4)  and 
the  hypotheses  of  Theorem  3.1  are  satisfied,  then 

(3.8)  $(sH*(s)  -  $*(s)i(((s)  t  0 

for  s  6  [o,5l  H  (5*, a*] .  ■ 

It  is  possible  to  relax  the  hypotheses  of  Theorem  3.1.  It  is  necessary 
only  that  z  and  z  be  piecewise  continuous. 

Theorem  3.3.  Suppose  that  (3.1)— (3.6)  are  satisfied  with  s  ♦  z(s),  z*(s) 

piecewise  continuous  and  having  points  of  discontinuity  i  -  l,...,n  and 

.  * 

i  -  l,...,n  ,  respectively.  If  there  exist  non-zero  constants  ki ,  i  » 
l,...,n  and  k*,  i  -  l,...,n*  such  that 


continuous  but ,  at 
solved  forward , 
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(3.9a)  4>(C^)  * 

(3.9b)  4>(C*)  =  k^C"), 

(3.9c)  x(C^)  3  kix(C1), 

for  i  =  l,...,n,  and  such  that 

(3.10a)  $*(5*~>  =  k*4>*(5*+), 

(3.10b)  /(£*')  =  k*/(q+) 

*  *—  v  *  *  X  *+  V 

(3.10c)  X  (q  )  =  klX  (c  ) 

for  i  -  1 , . . . ,n  ,  then  formula  (3.7)  continues  to  hold,  as  do  the 
conclusions  of  Corollary  3.2. 

*  *  *  * 


The  collections  of  functions  {($>,\{<,x.z} ,  C 4>  ,i|>  ,X  ,z  },  breakpoints 
*  * 

(C t } »  •  and  multipliers  {k^} ,  {k^}  satisfying  the  conditions  of 

Theorem  3.3  comprise  a  factorizations  of  the  two  point  boundary  value  problem 
(2.4). 

Theorem  3.3  is  central  to  the  viability  of  a  computer  implementation  of 
a  factorization  procedure  based  on  the  solution  of  (3.1 )— ( 3.6)  followed  by  the 
solution  of  (3.7)  at  selected  target  points  at  which  the  values  of  u  and 

A  A 

v  are  desired.  The  conditioning  functions  z,  z  ,  the  breakpoints 

and  the  multipliers  k^ ,  k*  are  chosen  adaptively  by  the  computer  during  the 

computation  in  order  to  stabilize  the  factorization. 


3.2.  Stable  Factorizations 


Definition  3.4.  A  factorization  for  (2.4)  is  bounded  if  there  exist 
constants  0  <  X  <  A  <  ®  such  that 

(3.11)  X2  <  $2( s)  +  ^2(s)  <  A2 

r  *  i  * 

for  all  s  €  1.0, a  1  with  a  similar  inequality  for  <j>  . 

Bounded  factorizations  are  stable  in  the  sense  of 

Theorem  3.5.  Let  $  X  anc*  ^  ^  »  X  be  the  solutions  of  the 

initial  value  problems  (3.1)- (3. 6)  under  the  relaxed  hypotheses  of  Theorem  3.3 
and  such  that  the  right  hand  sides  have  been  perturbed  by  functions  of 
magnitute  T  (in  the  norm).  Then  u,  v  solving  (3.7)  with  $ 

replacing  <p,  etc.,  are  solutions  of  a  perturbed  problem  (2.4),  for  which  the 
perturbations  of  the  input  data  are  of  magnitude  c(X,A)t  where  c(X,A)  is 
independent  of  the  problem.  a 

The  significance  of  Theorem  3.5  lies  in  the  fact  that  c(X,A)  is 
independent  of  the  data  of  the  problem  (e.g.  of  e  in  (2.1a)),  and  in  the 
fact  that  the  numerical  solution  of  the  initial  value  problems  (3. 1)— (3.6)  can 
be  interpreted  as  the  exact  solution  of  perturbed  equations.  The  magnitudes 
of  the  perturbations  of  the  equations  are  under  the  direct  control  of  the 
solution  tolerance,  for  example,  of  an  adaptive  solver.  The  adaptive 

it  * 

selection  of  the  functions  z,  z  the  breakpoints  5^,  ^  and  the 

multipliers  k^ ,  k£  by  the  computer  itself  is  essential  to  ensure  that  A/X 
and  A  are  of  reasonable  magnitude . 
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Example  3 .6 .  Normalized  Factorization 
Set 


(3.12) 


-(s)  ,  a<  s)tfr^(  s)+(  b(s)+c(s))<fr(s)i);(s)+d(s)i|)~(s) 

d>2(  s)  +  ^2(  s) 


It  is  not  hard  to  see  that  satisfying  (3.1)  and  (3.  with  this 

choice  of  z  have  the  property  that 


(3.13) 


>2(s)  +  ij)2(  s)  -  1 


2  2 

if  a  +3  «  1  (which  normalization  we  may  assume  without  loss  of  general¬ 
ity).  In  this  case  z  is  continuous  and  the  sets  of  breakpoints  and 
multipliers  are  empty. 

Example  3.7 .  Riccati  Factorization 

It  is  possible  to  normalize  the  boundary  conditions  (2.4b)  so  that 
either  a  *  1  and  ] 0 |  <  1  or  0*1  and  |a|  <  1.  If  a  =  1  set 

(3.14)  z(s)  =  a(sH(s)  +  b(s)\J;(s). 

In  this  case  it  follows  that  $  satisfying  (3.1)  is  constant, 

(3.15)  Hs)  =  l, 
and  equations  (3.2)  and  (3.3)  take  the  forms 

ij>'(s)  -  -c(s)  +  (a( s )-d( s) )i(i(  s)  +  b(s)'p"(s) 

X'(s)  -  f(s)  +  g(s)ij>(s)  +  a(s)x(s)  +  b(  s)^(s)x(s) . 


(3.16) 

*'(s) 

(3.17) 

X'(s) 

If  0*1  set 


I 


(3.18) 


z(s)  =  d(  s)i)/(  s)  +  c(s)<j>(s). 


In  this  case  it  follows  that  \|i  satisfying  (3.2)  is  constant. 


(3.19) 


ty(s)  =  1, 


and  equations  (3.1)  and  (3.3)  take  the  forms 


(3.20) 


(3.21) 


(j>'(s)  =  b(s)  +  (d(s)-a(s))<j>(s)  -t*  c(s)4>"(s) 


X'(s)  =  g(s)  +  f(  s)  <J)(  s )  +  d(  s)  x(  s)  +  c(  s).j)(s)x(s) . 


The  constraints  (3.15)  or  (3.19)  are  made  explicitly  in  (3.1),  (3.2),  and 
(3.3).  This  reduces  by  one  the  number  of  initial  value  problems  to  be  solved 
in  each  direction  and  the  lower  bound  \  equals  1.  Unfortunately  the 
solutions  of  the  resulting  (Riccati)  equations  may  not  exist  on  the  entire 
interval  [a, a  ].  But  if  they  do  not,  this  fact  is  signalled  by 

9  9  9  2  2  7 

$"(s)  +  \p  (s)  *  «  (i.e.,  1  +  ^“(s)  >  A  or  b  (s)  +  1  >  A")  in  which  case  a 
breakpoint  can  be  defined  and  the  solutions  renormalized.  That  is,  the 
factorization  algorithm  can  switch  adaptively  between  the  ^-problem  and  the 
ili-problem . 


i 


/ 


11 


4.  DIRECT  SOLUTION  OF  THE  MODEL  PROBLEM 

The  direct  method  of  solution  of  (2.1)  is  based  upon  the  discretization 
of  the  problem  by  a  finite  difference  (or  finite  element)  scheme  and  the 
subsequent  solution  of  the  resulting  matrix  equations  for  the  values  of  the 
discrete  solution.  The  dicretization  may  be  carried  out  in  the  second-order 
form  (2.1)  or  in  the  transformed  first-order  forms  (2.3)  or  (2.4).  In  either 
case,  a  complete  analysis  of  the  algorithm  should  consider  the  underlying 
numerical  processes  which  result  from  the  solution  of  the  discrete  matrix 
equations.  Let  us  carry  out  such  an  analysis  for  a  particular  discretization 
based  on  (2.4). 


4.1.  Block-Diagonal  Matrix  Equations 

A  difference  scheme  for  poblems  of  the  type  (2.4)  with  a(s)  »  d(s) 
has  been  proposed  in  [1],  where  such  problems  are  said  to  be  essentially 
diagonally  dominant.  In  the  spirit  of  the  particular  form  (2.3)  let  us 
restrict  ourselves  to  the  case  d( s)  =  0. 

Suppose  that  a  mesh 


o  = 


.00 


<  o 


(h) 


,(h) 


=  a 


with  it  c  has  been  given.  The 

of  the  partition;  e.g.,  h  =  max( 
proposed  in  [1]  can  be  written 


parameter  h  characterizes  the  fineness 
-0^h^).  Then  the  difference  scheme 


(h)  (h) 

Vl  ui 


KtiaAh)  +  civih))  +  +  will) 


v<h>-v<h) 
i+1  i 


- ± -  -  JL  (b  u(h)  +  b  U(h))  +  —  (g  +7  ) 

1  2  '  i+1  i+I  ii  '  2  V*i+1  V 


(4.1a) 


12 


where 


(4.2) 


0  if  |a1h1|  >  1,  at  <  0 


\  if  |aih1|  <;  1 


1  if  j  a± ^ |  >  1,  a±  >  0. 


Written  as  a  system  of  linear  equations  for  the  vector 

=*  [•••  u^  ,v^  ,ui+i,vi+i  ***J  n°dal  variables,  (4.1)  has  the  form 

(4.3)  A(h)U(h)  -  F(h) 

where  the  structures  of  and  are  given  in  Table  4.1,  based  on  the 

quantities 


--rV 


-i, 

^ib 

2  bi+l  ’ 


(4.4) 


t . 
l 


1, 

hi 

2  ^gi+gi+P’ 


Ci  "  "(1+Kihiai) * 


wi  “  ■''ihi0!’ 


Xi  -  1  "  (l-»c1)h 


1*1+1’ 


't  ”  “( 1-,ci^hici+i  * 


for  i  -  0,l,...,ra(h)-l, 
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Perform  forward  Gauss  elimination  with  pivoting  on  in  order  to 

bring  it  into  the  upper  triangular  form  shown  in  Table  4.2.  A  similar 
backward  elimination  beginning  with  the  last  row  of  A^b^  brings  it  into  the 
lower  triangular  form  shown  in  Table  4.3.  Observe  that  the  nodal  values  of 
the  solution  vector  then  satisfy  the  local  2x2  matrix  equations 

(4.5) 


* 

**i+l 


'i+lj 


- 

r  (h) 

X, 

u 

i 

v(h) 

* 

Lvi  J 

_xi+l 

for  i  ■  0, 1 , . . . ,ra^ .  Clearly  (4.5)  is  a  discrete  analog  of  (3.7).  However 

*  *  * 

now  the  entries  ^ and  X^  satisfy  explicit  one-step 

recursion  formulas  which  are  determined  by  the  nature  of  the  elimination 
algorithm. 

As  an  example ,  it  is  possible  to  perform  the  elimination  in  such  a  way 
that  one  or  the  other  of  the  following  two  recursions  hold : 


i|<  recursion: 


(4.6a)  $ 


i+1 


1, 


(4.6b)  ij; 


i+1 


•■J^d+Kjhja^+KjhjCj-U-  ~  b14/1)(l-K1)h1c1 

h|  iT 

"  T  bi+i^i(1+Kihiai)“’cihici)+(1"  T  Vi)(1-<1_'<i)hiai+1) 


1 

~  (gQ'^1+b1X1)('l'1(l+'<lh1al)-«lhlc1)+x1(l+Kihia1) 
i+1  “  ~~h~  hT  5 

‘  T  W^1+Kihiai)',cihici)+(1"  y  b^iXi-n-^JVi+i) 


(4.6c)  x 
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ifr  recursion; 


(4. 7a)  *i+1 


bi+1(1+Kihiar<thicl*i)+(*r  t  b1)(l-»-«i)hi*l4.i) 

hi 

(l+ic^apc^c^H^-  —  b1)(l-<i)hici+1 


(4.7b)  *±+1  -  1, 


(4.7c) 


4+1 


i 

(gi^i+l)+xi)+xi(1+^Lhiai_Kihicil,i)+((,’i"  T  bi^ihlcixi 

(l+Kihiai-<1hic14>1)-(<(li-  j-  bi](l-<i)hic1+1 


It  Is  readily  seen  (by  letting  h  0)  that  the  explicit  recursion 
formulas  (4.6)  and  (4.7)  can  be  regarded  as  special  one-step  methods  for  the 
initial  value  problems  (3. 15)— (3— 17)  and  (3.19 )— (3— 21 ) ,  respectively.  That 
is,  the  closure  of  (4.6)  is  (3. 15)— (3— 17)  and  the  closure  of  (4.7)  is  (3.19)- 
(3-21).  Of  course,  it  is  also  possible  to  devise  an  elimination  scheme  for 
(4.3)  which  yields  the  normalized  factorization  (z  given  by(3.12))as  its 
closure.  Analogous  recursions  can  be  given  for  the  backward  elimination. 
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5.  STABLE  NUMERICAL  PROCESSES 

The  explicit  one-step  recursion  formulas  (4.6)  and  (4.7)  are  examples  of 
a  general  one-step  numerical  process  (cf.  [2])  of  the  form 

(5.1)  ^i+1  ™  $^(y^,h^),  ^  *  0,1,2,... 

where  y^.G^  €  R‘  .  We  have  already  noted  that  (4.6)  and  (4.7)  are.  consistent 
with  (i.e.,  the  closures  of)  particular  realizations  of  the  initial  value 
problems  (3.1 )— ( 3.3).  They  are  also  stable. 

Definition  5.1  (Stability).  The  process  (5.1)  is  stable  if  it  is  bounded, 
(5.2a)  3  y1 B  <  K, 

and  if  there  exist  a  neighborhood  B{(y)  *  (z:  8z~y  ! q<  6}  of  yQ,  and  an 

L  <  1  such  that 

(5.2b)  ,zi+l-yi+ll  <L|zi_ui11 


for  given  by 


zi+1  -  wv 


with  zQ  6  B6(yQ). 

Definition  5.1  is  an  adaptation  of  BN-statibility  (cf.  [4]).  In 
contrast  to  the  concepts  of  A-stability  and  B-stability,  the  definition  is 
made  without  reference  to  a  particular  test  equation.  If  L  <  1  the  process 
is  strongly  stable. 

Consider  the  model  problem  (2.1)  with 

a(  s)h 
e 


(5.3a) 


<  -1 
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and 


(5.3b) 


b(s)  >  0. 


Under  these  conditions  it  is  not  hard  to  see  that  recursion  (4.6)  will  hold 
with  -  0  in  (4.2).  In  fact,  (4.6b)  reduces  to 


hibi 


(5.4) 


Ki+i 


,2, 

h.b.a.  , 
r  i  l  i+1 


‘WW 


)h  +  (e-hiai+i) 


which  can  be  shown  to  be  an  order  1  method  for  (3.16). 

Theorem  5.2.  For  c  <  e^,  with  eQ  such  that  conditions  (5.3)  hold,  the 
recursion  (5.4)  remains  in  the  semi-infinite  strip 


(5.5) 


2  - 


4>  < 


r^-2.)  “ 

>■  ah  J 

“TB — ’ 


a  >  maxja(s)|,  b  >  max|b(s)j 


and  is  strongly  stable,  independently  of  e  .  ■ 

Results  analogous  to  Theorem  5.2  hold  for  the  x  recursion  (4.6c),  and 

1 h3( S ) I 

also  for  the  recursions  (4.7)  in  the  event  that  | — jr-4|  >  1.  For  the  case 
i hs( s ) i 

| - ^-|  <1  it  is  not  possible  to  say  a  priori  whether  recursion  (4.6)  or 

(4.7)  will  be  used.  It  is  likely  that  the  matrix  reduction  algorithm  will 
switch  between  them,  depending  upon  the  behavior  of  a(s),  b(s),  and  on  the 

pivoting  strategy  in  a  particular  case.  It  is  not  our  intent  to  give  an 
exhaustive  analysis  here. 

As  a  final  example  to  illustrate  the  ideas  discussed  above,  note,  that 
another  first  order  accurate,  stable,  explicit  recursion  for  the  initial  value 


20 


problem  (3.2)  can  be  obtained  by  applying  the  implicit  Euler  method, 


(5*6)  *i+i  "  *i  +  hi  I"  z  +  (“T^i+i  +  bi+i*i+iJ* 


whence 


e_hiai+l  *ehibi+1*®*i"h1)  l/2 


Ki+1 


2ehibi+l 


{1  - 


( e~hi ai+l ) * 


} 


Thus 


(5.7) 
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6.  CONCLUSION 

We  have  shown  on  a  model  problem 

1)  that  there  exist  factorizations  of  the  original  two  point  boundary 
value  problem  into  systems  of  initial  value  problems  with  guaranteed 
stability  properties; 

2)  that  there  can  be  many  methods — some  specially  designed — which  solve 
these  initial  value  problems.  Some  of  these  methods  will  produce 
the  identical  numerical  processes  as  the  finite  difference  (or 
finite  element)  methods  applied  to  the  original  boundary  value  value 
problems . 

In  view  of  these  two  points,  we  suggest  that  it  may  be  advantageous  to 
avoid  the  finite  element  or  finite  difference  approach  altogether,  and  rather 
to  study  the  class  of  stable  factorizations  directly,  with  the  goal  of 
selecting  both  an  optimal  factorization  for  the  given  boundary  value  problem, 
and  also  an  optimal  method  (possibly  a  special  one)  for  solving  the  initial 
value  problems  of  the  factorization.  In  this  way  one  might  employ  the  modern 
ideas  and  adaptive  approaches  widely  used  in  today's  software  for  the 
numerical  treatment  of  initial  value  problems  for  ordinary  differential 
equations . 


» 

i 
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